% Working m-file for model.mod

clear
close all

%% Parameters
% Target values
%--------------------------------------------------------------------
hU_SS = 1;                      % average hours worked in SS
LU_SS = 20;                     % leverage in SS
PU_SS = 0.05/4;                 % crisis probability

% Households
%--------------------------------------------------------------------
beta_p = 1.02^(-.25);           % preference discount factor
nu_p = 1e+10;                       % Frisch elasticity of labor supply                  

% Firms
%--------------------------------------------------------------------
alpha_p   = 0.36;               % capital-income share
delta_p   = 0.05;              % depreciation rate

% Fund managers
%------------------------------------------------------------------
lam_p = 0.5;                    % fire sale price (no fire sale if lam_p = 1)

% Banks
%--------------------------------------------------------------------
chi1_p = 0.95;                  % Surviving probability
chi0_p = 0.05;                  % Dividend payout ratio

% TFP shock
%--------------------------------------------------------------------
rhoa_p = 0.9;
siga_p = 2/100;

%% Steady state
hU = hU_SS;  % hours worked
LU = LU_SS;  % leverage
PU = PU_SS;  % crisis probability
aU = 0;      % TFP shock

RbU = 1/beta_p;
ee = 1;
while ee > 1e-10
    RbUp = find_RbU(RbU,LU,PU,nu_p,alpha_p,siga_p,lam_p,delta_p,beta_p);
    ee = abs(RbU-RbUp);
    RbU = RbUp;
end
[~,gam_p,rkhatsU,rkhatU] = find_RbU(RbU,LU,PU,nu_p,alpha_p,siga_p,lam_p,delta_p,beta_p);

RkU = exp(rkhatU)+1-delta_p;
kU = (alpha_p/(RkU-1+delta_p))^(1/(1-alpha_p));
yU = kU^(alpha_p);
iU = delta_p*kU;
cU = yU - iU;
psi_p = (1-alpha_p)*yU;
nU = kU/LU;
n0_p = (1-chi1_p*(chi0_p*((alpha_p*yU/kU+1-delta_p)*LU-RbU*(LU-1)-1)+1)).*nU;
xiU = exp(rkhatU)/exp(rkhatsU);

%% Check equations
murkU = log(alpha_p*((1-alpha_p)/psi_p)^((1-alpha_p)/(alpha_p+1/nu_p)))+(1+nu_p)/(alpha_p*nu_p+1)*rhoa_p*aU-(1-alpha_p)/(alpha_p*nu_p+1)*log(kU);
sigrk_p = (1+nu_p)/(alpha_p*nu_p+1)*siga_p;
PU = 1/2.*(1+erf((rkhatsU-murkU)/(sigrk_p*sqrt(2))));
pi = 3.141592653589793;
ERkp1d = exp(murkU+sigrk_p^2/2)*1/2*(1+erf((rkhatsU-murkU-sigrk_p^2)/(sigrk_p*sqrt(2))))+(1-delta_p)*PU;
    
ee(1) = -1/(beta_p*RbU) + ((cU-psi_p*hU^(1+1/nu_p)/(1+1/nu_p))/(cU-psi_p*hU^(1+1/nu_p)/(1+1/nu_p)))*(ERkp1d*LU/(RbU*(LU-1))-lam_p*PU+1-PU);

ERkp1u = exp(murkU+sigrk_p^2/2)*(1/2)*(1+erf((murkU+sigrk_p^2-rkhatsU)/(sigrk_p*sqrt(2))))+(1-delta_p)*(1-PU);
ee(2) = -RbU*(1-PU) + ERkp1u - lam_p*RbU^(2)*(1-gam_p)*(1+lam_p*(1-gam_p))*(LU-1)/(LU^2)/(RbU*(1-1/LU)*(1+lam_p*(1-gam_p))-1+delta_p)*1/(sqrt(2*pi)*sigrk_p)*exp(-1/2*((rkhatsU-murkU)/sigrk_p)^2);
